home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / DE118I.ZIP / RELTIV.C < prev    next >
C/C++ Source or Header  |  1991-11-04  |  3KB  |  131 lines
  1. #ifndef NOINCS
  2. #include "mconf.h"
  3. #include "prec.h"
  4. #include "const.h"
  5. #include "ssystem.h"
  6. #endif
  7.  
  8. /* Relativistic corrections.  Reference:
  9.  *
  10.  * Newhall, XX, E. M. Standish, and J. G. Williams, "DE102: a
  11.  * numerically integrated ephemeris of the Moon and planets
  12.  * spanning forty-four centuries," Astronomy and Astrophysics
  13.  * 125, 150-167 (1983).
  14.  */
  15.     
  16. extern DOUBLE C;
  17. extern DOUBLE GMs[];
  18. extern DOUBLE Rij[NTOTAL][NTOTAL];
  19. extern DOUBLE Rij3[NTOTAL][NTOTAL];
  20. extern DOUBLE vnewt[6*NTOTAL];
  21. DOUBLE v2j[NTOTAL];
  22.  
  23. reltiv( y, v )
  24. DOUBLE y[], v[];
  25. {
  26. DOUBLE xd, yd, zd, xs, ys, zs, xv, yv, zv, csqi;
  27. DOUBLE rc, r, s, temp, rci;
  28.  
  29. int i, j, k, ii, jj;
  30.  
  31. csqi = One/(C*C);
  32.  
  33. /* velocity terms */
  34. jj = 6*FMASS;
  35. for(j=FMASS; j<NTOTAL; j++)
  36.     {
  37.     xv = y[jj];
  38.     yv = y[jj+2];
  39.     zv = y[jj+4];
  40.     v2j[j] = (xv * xv + yv * yv + zv * zv) * csqi;
  41.     jj += 6;
  42.     }
  43.  
  44. ii = 6*FMASS;
  45. for(i=FMASS; i<IAROIDS; i++)
  46.     {
  47.     v[ii] = Zero;
  48.     v[ii+2] = Zero;
  49.     v[ii+4] = Zero;
  50.     rci = Zero;
  51.     for(k=FMASS; k<NTOTAL; k++)
  52.         {
  53.         if( k == i )
  54.             continue;
  55.         rci += GMs[k]*Rij[i][k];
  56.         }
  57.     xs = Zero;
  58.     ys = Zero;
  59.     zs = Zero;
  60.     jj = 6*FMASS;
  61.     for(j=FMASS; j<NTOTAL; j++)
  62.         {
  63.         if( j == i )
  64.             {
  65.             jj += 6;
  66.             continue; /* skip to next j if i = j */
  67.             }
  68.  
  69.         xd = y[jj+1] - y[ii+1];
  70.         yd = y[jj+3] - y[ii+3];
  71.         zd = y[jj+5] - y[ii+5];
  72.  
  73.         rc = -Four * csqi * rci;
  74.  
  75.         s = Zero;
  76.         for(k=FMASS; k<NTOTAL; k++)
  77.             {
  78.             if( k == j )
  79.                 continue;
  80.             s += GMs[k]*Rij[j][k];
  81.             }
  82.         rc -= csqi * s;
  83.  
  84.         rc += v2j[i];
  85.         rc += Two * v2j[j];
  86.  
  87.         r = y[ii] * y[jj] + y[ii+2] * y[jj+2] + y[ii+4] * y[jj+4];
  88.         rc -= Four * csqi * r;
  89.  
  90.         s = -xd * y[jj] - yd * y[jj+2] - zd * y[jj+4];
  91.         s *= Rij[i][j];
  92.         rc -= OneandaHalf * csqi * s * s;
  93.  
  94.         s = xd * vnewt[jj] + yd * vnewt[jj+2] + zd * vnewt[jj+4];
  95.         rc += Half * csqi * s;
  96.  
  97.         temp = GMs[j] * Rij3[i][j];
  98.  
  99. /*
  100.  *        rc = temp * (One + rc );
  101.  */
  102.         rc = temp * rc;
  103.         xs += xd * rc;
  104.         ys += yd * rc;
  105.         zs += zd * rc;
  106.  
  107.         s = -xd * (Four*y[ii] - Three*y[jj])
  108.             - yd * (Four*y[ii+2] - Three*y[jj+2])
  109.             - zd * (Four*y[ii+4] - Three*y[jj+4]);
  110.         s = s * temp * csqi;
  111.         xs += s * (y[ii] - y[jj]);
  112.         ys += s * (y[ii+2] - y[jj+2]);
  113.         zs += s * (y[ii+4] - y[jj+4]);
  114.  
  115.         temp = ThreeandaHalf * csqi * GMs[j] * Rij[i][j];
  116.         xs += temp * vnewt[jj];
  117.         ys += temp * vnewt[jj+2];
  118.         zs += temp * vnewt[jj+4];
  119.         jj += 6;
  120.         }
  121.     v[ii] = xs;
  122.     v[ii+2] = ys;
  123.     v[ii+4] = zs;
  124.     v[ii+1] = y[ii];
  125.     v[ii+3] = y[ii+2];
  126.     v[ii+5] = y[ii+4];
  127.     ii += 6;
  128.     }
  129. }
  130.  
  131.